## This file assesses the matches with whites using only those who were shown their DA
## but allowing matches across different DAs within the design

options(digits = 4, width = 132)

library(here)
library(tidyverse)
library(RItools)

source(here::here("Design", "nonbimatchingfunctions.R"))

load(here::here("Design", "matches_DA_new.rda"), verbose = TRUE)
load(here::here("Data", "wrkdat_DA0_new.rda"), verbose = TRUE)

wdat0 <- wrkdat_DA0_new %>%
  dplyr::select(-geometry) %>%
  filter(!is.na(vm_change) & !is.na(social.capital01) & !is.na(community.resp01) & !is.na(vm)) %>%
  droplevels()
dim(wdat0)

wdat0 <- left_join(wdat0, matches_DA_new)
table(is.na(wdat0$pair))

## Easier to think about "higher perceiver" versus "lower perceiver" when it comes to asking whether the matching worked.
## Further simplify the working file
wdat1 <- wdat0 %>%
  filter(!is.na(pair)) %>%
  group_by(pair) %>%
  mutate(perc_more = rank(vm) - 1) %>%
  ungroup() %>%
  mutate(pairF = factor(pair)) %>%
  select(
    perc_more, da_prop_vm_20pct_06, vm_change, pair, pairF, social.capital01, community.resp01,
    vm, vcid, dauid, csd_pop_06, csd_prop_vm_20pct_06,
    da_pop_dens_06, csd_urban_06, csd_pop_dens_06, community_area_km, da_pop_06
  ) %>%
  droplevels()

xb1_ranked0 <- xBalance(perc_more ~ da_prop_vm_20pct_06 + vm_change,
  strata = list(unstr = NULL, pair = ~pair), report = "all", data = wdat1
)
xb1_ranked0$overall
xb1_ranked0$results

xb1_ranked <- balanceTest(perc_more ~ da_prop_vm_20pct_06 + vm_change + strata(pair), data = wdat1, p.adjust = "none")
xb1_ranked$overall
xb1_ranked$results[, c("Control", "Treatment", "adj.diff", "std.diff", "p"), ] # , "pair"]


## An alternative for the overall test which should be fairly close to it
library(formula.tools)
library(coin)
coin_fmla <- da_prop_vm_20pct_06 + vm_change ~  vm | pairF
coin_test <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic")
coin_test_perm <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic", distribution = approximate(nresample = 1000))
coin::pvalue(coin_test)
coin::pvalue(coin_test_perm)

coin_fmla_rank <- da_prop_vm_20pct_06 + vm_change ~ perc_more | pairF
coin_test_rank <- independence_test(coin_fmla_rank, data = wdat1, teststat = "quadratic")
coin_test_rank_perm <- independence_test(coin_fmla_rank,
  data = wdat1,
  teststat = "quadratic", distribution = approximate(nresample = 1000)
)
coin::pvalue(coin_test_rank)
coin::pvalue(coin_test_rank_perm)

pairdiffs <- wdat1 %>%
  group_by(pair) %>%
  summarize(
    da_prop_vm_20pct_06 = abs(diff(da_prop_vm_20pct_06)),
    vm_change = abs(diff(vm_change)),
    vm = abs(diff(vm)),
    da_pop_06 = abs(diff(da_pop_06)),
    csd_pop_06 = abs(diff(csd_pop_06)),
    csd_pop_dens_06 = abs(diff(csd_pop_dens_06)),
    csd_prop_vm_20pct_06 = abs(diff(csd_prop_vm_20pct_06)),
    # da_pop_dens_06 = abs(diff(da_pop_dens_06)),
    community.area.km = abs(diff(community_area_km))
  )

sapply(pairdiffs, summary)

sapply(pairdiffs, function(x) {
  quantile(x, sort(unique(c(seq(0, 1, .1), seq(.9, 1, .01)))), na.rm = TRUE)
})
#       pair da_prop_vm_20pct_06 vm_change     vm da_pop_06 csd_pop_06 csd_pop_dens_06 csd_prop_vm_20pct_06 community.area.km
# 0%     1.0          0.00000000  0.000000 0.0100       0.0          0            0.00             0.000000           0.01243
# 10%   23.3          0.00000000  0.001581 0.0300      23.1          0            0.00             0.000000           1.33334
# 20%   45.6          0.00009896  0.005265 0.0660      46.0      11500           24.10             0.007092           4.90736
# 30%   67.9          0.00044887  0.008501 0.1000      84.9      55164           56.39             0.014891           8.51434
# 40%   90.2          0.00091541  0.011661 0.1300     128.4     101246          166.82             0.031655          22.78940
# 50%  112.5          0.00143456  0.014479 0.1700     182.0     257821          282.19             0.043136          46.11770
# 60%  134.8          0.00243473  0.017898 0.2100     249.2     413837          399.00             0.072808          98.16405
# 70%  157.1          0.00357245  0.019897 0.2710     375.5     685513          858.66             0.130674         178.32007
# 80%  179.4          0.00541113  0.022633 0.3400     586.4    1142080         2074.33             0.191346         383.60896
# 90%  201.7          0.00901861  0.026221 0.4800    1122.8    1925240         2904.36             0.302753         862.84720
# 91%  203.9          0.00908558  0.026390 0.4993    1324.0    1925240         2912.35             0.306581        1177.05877
# 92%  206.2          0.00941058  0.026802 0.5100    1454.0    2016255         2931.77             0.328060        1214.61644
# 93%  208.4          0.00991629  0.027148 0.5278    1565.4    2124912         3090.47             0.331926        1530.49800
# 94%  210.6          0.01078593  0.027474 0.5824    1804.7    2235158         3118.48             0.336357        1739.70889
# 95%  212.9          0.01112658  0.028284 0.6740    2072.7    2300202         3193.78             0.369550        1948.84696
# 96%  215.1          0.01123025  0.028435 0.7632    2452.2    2390318         3303.52             0.377446        2048.72308
# 97%  217.3          0.01163234  0.029130 0.9117    2638.6    2428464         3393.22             0.399906        2658.12489
# 98%  219.5          0.01268062  0.029422 1.0362    3827.6    2444732         3464.38             0.416627        3909.28556
# 99%  221.8          0.01409258  0.029688 1.3231    6011.4    2479360         3902.14             0.433720        4998.60722
# 100% 224.0          0.01626694  0.029860 4.1900   11429.0    2499379         4537.65             0.499291       29227.81244

### facts for the body.tex

# number of people in final design
nrow(wdat1)
# [1] 448
# number of DAs
length(unique(wdat1$dauid))
# [1] 443

system("touch Design/match_assess_DA_new.done")
